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We examine one- and two-dimensional (ID and 2D) models of linearly coupled lattices of the 
discrete-nonlinear- Schrodinger type. Analyzing ground states of the systems with equal powers in 
the two components, we find a symmetry-breaking phenomenon beyond a critical value of the squared 
/^-norm. Asymmetric states, with unequal powers in their components, emerge through a subcritical 
pitchfork bifurcation, which, for very weakly coupled lattices, changes into a supercritical one. 
We identify the stability of various solution branches. Dynamical manifestations of the symmetry 
breaking are studied by simulating the evolution of the unstable branches. The results present 
the first example of spontaneous symmetry breaking in 2D lattice solitons. This feature has no 
counterpart in the continuum limit, because of the collapse instability in the latter case. 



I. INTRODUCTION 

Dynamical lattices and their applications have become 
an area of increasing interest over the past decade, as 
shown by a multitude of recent reviews of the topic [1]- 
[2,]. This growth was driven by a wide array of physical 
realizations, in fields as diverse as light propagation in 
optical waveguide arrays dynamics of Bose-Einstein 
condensates (BECs) in periodic potentials (optical lat- 
tices) [5], micro- mechanical cantilever arrays models 
of DNA |6|], and others. A key model that has been 
widely used and analyzed in each of the above areas is 
the discrete nonlinear Schrodinger (DNLS) equation [2j]. 
In these applications, it emerges either as a tight-binding 
approximation (as in the case of BECs trapped in optical 
lattices), or at the level of an envelope- wave expansion of 
the underlying physical field (such as the electromagnetic 
field of light in the optical systems). 

One aspect of this class of discrete dynamical mod- 
els which remains perhaps less explored concerns their 
multi-component generalizations, which are relevant to 
many fields where dynamical lattices are natural mod- 
els. For instance, in the case of waveguide arrays, one 
may consider settings with two orthogonal polarizations 
of light [7], or two different wavelengths, see e.g., Ref. [8]. 
Similarly, in the BEC context, one may consider multi- 
species condensates in the form of mixtures of different 
hyperfine states in ^^Rb [i, [T^ and ^^Na [11], or mix- 
tures of different atomic species, such as Na-Rb, K-Rb, 
Cs-Rb, Li-Rb, as well as Li-Cs (see, e.g., Ref. |12|] and 
references therein). 

While the above settings are typically modeled by sys- 
tems of DNLS equations which are coupled by nonlinear 



terms, such as the ones accounting for the cross-phase 
modulation in optics, or collisions between atoms be- 
longing to different BEC species, it is also relevant to 
consider linearly coupled DNLS equations (or, in other 
discrete settings, linearly coupled Ablowitz-Ladik equa- 
tions [HI). In the optics context, such systems of linearly 
coupled DNLS equations are relevant to various applica- 
tions: for example, linear coupling occurs among the two 
modes inside each waveguide, which may be induced by 
a twist of the waveguide (for linear polarizations), or by 
birefringence (for circular polarizations), or in a dual-core 
structure of the waveguide [8]. On the other hand, in 
BECs, a linear coupling may be imposed by an external 
microwave or radio-frequency field, which can drive Rabi 
[T3 | or Josephson [l5| oscillations between populations of 
two different states. 

In the present work, motivated by the earlier works 
m, [13, Us 9 which studied effects of the linear coupling 
in various continuum optical models, and also the analy- 
sis of the symmetry breaking of discrete two-component 
solitons in the linearly coupled Ablowitz-Ladik system 
P^], we consider a system of two DNLS equations cou- 
pled solely by linear terms. In the optical setting, this 
would be the above-mentioned discrete analog of a dual- 
core fiber [16]. Such a model has also been proposed 
as a means for realization of all-optical switching in ul- 
trashort photonic-crystal couplers [T9]. In BECs, it may 
model the dynamics of two coupled hyperfine states (e.g., 
of ^^Rb), where we imply the use of Feshbach-resonance 
techniques [20] to nullify the nonlinear interaction be- 
tween the components (by rendering the respective inter- 
species scattering length equal to zero). At the same 
time, two spin states may be linearly coupled due to a 
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FIG. 1: Bifurcation diagram of solutions for the anti- 
continuum limit, e = 0; r and E are the asymmetry parameter 
and the half of the total squared norm, respectively (see the 
text for the mathematical definitions). The solid and dashed 
lines show stable and unstable solutions, respectively. 



0.5 



1.55 



1.5 



1.45 



1.2 



1.3 

1^ 



1.4 



1.5 



-0.5 



3.5 



4.5 



FIG. 2: Bifurcation diagrams for e = 1.6 in the ID model. 
The dashed-dotted line indicates solutions found by means of 
the variational approximation, while solid and dashed lines 
show, respectively, numerically found stable and unstable 
steady- state solutions. 



resonant spin-flipping radio- frequency field, as mentioned 
above. In both media (optical and atomic), the relevant 
model takes the following form: 



iUt = KeA2U^KV^\UfU 
iVt = KeA2V ^KU ^ \V\^ V 



(1) 



where U = U{x^t) and V = V{x^t) are the wave func- 
tions in BEG or electric field envelopes in optics {x is 
realized as a set of discrete coordinates), A2 is the dis- 
crete Laplacian, formed by the centered difference in each 
of the relevant dimensions, K is the strength of the lin- 
ear coupling between fields U and and e determines 
the couplings between adjacent sites of the lattice. Note 
that, for convenience, the full coupling constant is de- 
fined as Ke; this convention will allow us to eliminate K 
from the analysis presented below. 

Our aim is to find and explore in detail the symmetry- 
breaking bifurcation of the ground-state single-pulse so- 



FIG. 3: The top panel shows the bifurcation diagram in the 
2D model for e = 0.25, in the same way as the ID diagram is 
shown in Fig. [2] The bottom panel displays the dependence of 
the solution's squared norm, upon the chemical potential, 
/i, for the symmetric solutions. Unlike the ID case, there are 
two different symmetric solutions for many values of E, re- 
sulting in both stable and unstable solutions for norms below 
the value at which the symmetric and asymmetric solution 
branches intersect. 



lution of Eqs. ([T]), similar to the earlier studies performed 
for continuum models in Refs. |lB, [13, We will ad- 
dress this problem analytically, by means of a variational 
approximation (VA), and numerically, via computations 
of the steady-state bifurcations and linear stability, as 
well as through direct numerical simulations testing the 
stability or instability of states under consideration. In 
this way, we obtain a complete bifurcation diagram of the 
discrete model. A unique advantage of performing this 
analysis in the discrete setting is that the bifurcation dia- 
gram can be produced not only for one-dimensional (ID), 
but also for two-dimensional (2D) lattices. The latter is 
impossible in usual continuum models of the cubic-NLS 
type because of the collapse instability [21]. However, 
it was demonstrated in Ref. ^ (see also Ref. ^) that 
the DNLS with sufficiently weak inter-site coupling [i.e., 
small coefficient e in Eqs. ([T|)] gives rise to stable discrete 
2D solitons. This permits us to construct a bifurcation 
diagram for the 2D lattice and compare it to the ID 
counterpart. To the best of our knowledge, this is the 
first time that a symmetry-breaking bifurcation is found 
and analyzed in a 2D model with the cubic nonlinearity. 
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FIG. 4: Plots of solutions from the branch in Fig. [2] for E — 3.4. The top-row figures show the solution profiles found by 
means of the numerical (Un, Vn) and variational ("analytical", Ua, Va) methods. The bottom row plots illustrate the stability 
eigenvalues for the numerical solution. The first column presents a stable stationary asymmetric solution from the outer (upper) 
branch in Fig. [21 the second column is an unstable asymmetric solution of the inner branch, and the last column is taken from 
the stable part of the family of symmetric solutions, with r = 0. 



The paper is structured as follows. In Section 2, we 
present the model and analytical results, based on the 
variational method. In Section 3, numerical results are 
reported. Finally, in Section 4, we summarize our find- 
ings and discuss directions for future work. 



II. THE MODEL AND ANALYTICAL 
CONSIDERATIONS 

In what follows, we seek steady- state (standing- wave) 
solutions to system ([1]). We use the standing- wave 
ansatz. 



U (f , t) = VK u (f ) exp[-iK (/i - 2De) t] 

V{x,t) = VKv{x) exp[-iK{/j.-2De)t] 

where u{x) and v{x) are real-valued functions, /i is the 
chemical potential (propagation constant) in the BEG 
(optics) model, shifted by the constant D; the latter takes 
the values D = 1 and 2 for the ID and 2D cases, respec- 
tively. Substituting these expressions in Eqs. ^ leads to 
stationary equations, which, in the ID model, take the 
form: 



lilUn = eAiUn + + 



(2) 



Where AiWn = i^n+i + ujn-i- In the 2D case, the sta- 
tionary equations are: 
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flVn 



(3) 



Where A2^n,m ^n+l,m + Wn-l,m + ^n,m+l + ^n,m-l. 

Below, we aim to construct symmetric {u = v) and 
asymmetric {u ^ v) solutions of Eqs. ([2|) and ([3|). In 
order to develop an analytical approximation to the so- 
lutions, we resort to the variational method [23]. To this 
end, we notice that Eqs. ([2]) and (|3]) can be derived from 
the following Lagrangians: 
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FIG. 5: Time-evolution plots for perturbed one-dimensional solutions from Fig. [2] The top (middle) row depicts the space-time 
contour plot evolution of \U\'^ (1^1^)- The first two columns are for two different perturbations of the solution shown in the 
second column of Fig. |4l the first column is generated by a perturbation which pushes the solution towards the upper stable 
asymmetric solution branch of the bifurcation diagram, while, in the second column, the perturbation pushes the solution 
towards the r = solution. The last column is for a perturbation of the unstable symmetric solution (on the r = curve) at 
E = 4.1, which pushes the solution towards the upper stable solution branch. The bottom row of figures displays the value of 
r as the solutions evolve in time, shown by the solid line, and the constant value of r for the steady states at the same energy 
level (dashed horizontal lines). 



Further, we employ natural ansdtze^ {un^^n} = 
{A,B}e-^\^\ and K,m,^n,m} = {A, 5}e-^l-le-^l-l, 
with free constants A, B and A > 0, in the ID and 
2D cases, respectively. This choice is motivated both 
by the exponential decay of the solutions' tails far from 
the soliton's center and by the fact that only this type of 
the trial functions makes the approximation straightfor- 
wardly tractable [24]. Plugging the ansdtze into Eqs. (j3|) 
and ([5]) and analytically evaluating the resulting geomet- 
ric series, we arrive at the following expressions for the 



effective Lagrangians: 



^ID 



L 



2D 



AB 



AB 



coth A 



B"*) coth (2A) + €{A^ + B^) cosechA (6) 



B' 



coth^ A 



B'^ ) coth^ 2 A 



2e {A^ + (cosechA) coth A. 



(7) 
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FIG. 6: Cross-section plots of the asymmetric solutions belonging to corresponding curve in Fig.[3l for E = 1.435. The top- row 
figures show the solutions found by means of the numerical (Un, Vn) and variational ("analytical", Ua, Va) methods, and the 
bottom row plots illustrate the stability eigenvalues for the numerical solution. As can be seen, the first and second columns 
represent, respectively, stable and unstable solutions belonging to the upper (outer) and inner curves (asymmetric branches) 
of the bifurcation diagram, respectively. 



The static form of the Euler-Lagrange equations follow- 
ing from here, 9i^iD,2D/S (A, A, B) = 0, is 

I {A^ + B^) cosech^A - ^ (^^ + cosech22A 

-A^cosech^Ae [A^ + B'^) cosechAcoth A = 0, 

-/iAcothA + A^coth2A + 5coth A + 2eAcosechA = 0, 

-/i5cothA + 5^coth2A + AcothA + 2e5cosechA = 0, 
for the ID case, and 

/i (A^ + B'^) coth Acosech^A 

- {A^ + B^) coth 2Acosech^2A - 2AB coth Acosech^ A 

- 2e {A^ + B^) (cosechAcoth^ A + cosech^A) = 0, 

- /I A coth^ A + coth^ 2 A + 5 coth^ A + 
4e74cosechA coth A = 0, 

- liiB coth^ A + 5^ coth^ 2A + A coth^ A 

+4e5cosechA coth A = 0, 

for the 2D case. These equations were solved for A, 5, 
and A and will be compared, in the next section, to the 
full numerical solution of Eqs. ([2]) and (j3|). 



Another analytically tractable case corresponds to the 
ant i- continuum limit of e = 0. For the symmetric branch, 
we then have Un = Vn = or Un = Vn = V/i — 1, while for 
the asymmetric branch, one needs to solve a system of al- 
gebraic equations, jJ.Un = Vn^U^, J^^l = U^^UnVn^V^^. 

The solution is shown in Fig.[Tl which displays the respec- 
tive symmetry-breaking bifurcation by means of a plot 
of the asymmetry measure, r = {Ei — E2)/{Ei + £^2), 
versus half the total norm, E = {Ei + £^2)/2, where 
{El, E2} = Xln {^n^^n} uorms of the two compo- 

nents of the solution (in the BEC model, they are propor- 
tional to the number of atoms in the two atomic states, 
while in the optical setting they measure the total power 
of the beams in the two coupled lattices). As seen from 
the figure, the pitchfork bifurcation is supercritical in this 
limit (this will be compared to typical behavior for finite 
e below). 



III. NUMERICAL METHODS AND RESULTS 

Numerical solutions to Eqs. ([2]) and (|3]) were found by 
usi ng t he method of the pseudo-arclength continuation 
[25I . [2^ . Another typical approach to solving systems of 
nonlinear equations for various values of a control param- 
eter relies upon parameter continuation; however, this 
method fails for solution-parameter pairings where the 
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FIG. 7: Same as Fig. (|6)) for two symmetric solutions found at ^ = 1.435. 



resulting Jacobian is singular. On the other hand, the 
pseudo-arclength continuation addresses this problem by 
introducing an extra pseudo-arclength parameter, and in- 
cluding an additional equation into the system, which 
makes the solution and control parameter dependent 
upon the pseudo-arclength. Calling the pseudo-arclength 
parameter 5, the additional equation, F{u^v^ ji^ s) = 0, 
must be chosen such that F{u^v,fl^s = 0) = 0, where 
{u^v^p) is a solution of Eqs. ([2]) or (|3]). Thus, for F we 
used F{u^ /i, s) = \u — -\- \v — + |/i — flf — s^. 

Once the steady states were identified, their stability 
was examined in the framework of linear stability analysis 
by substituting a perturbed solution, in the form of 



U{x, t) = e-^^^ [u{x) + a{x) e^^ + . 
V{x,t) = e-^^^ [v{x) + c{x) e^^ + d*{x) e^*^] , 



(8) 



in Eqs. ([T]), the asterisk standing for complex conjuga- 
tion. The linearized equations for the perturbation eigen- 
modes a, 6, c, d were solved numerically, yielding eigenval- 
ues A associated with them. 

The results are summarized in Figs. [2] and [3l for the 
ID and 2D cases, respectively. In the ID case, the sym- 
metric branch, with r = 0, is stable for E < 3.82. Be- 
yond this critical point, it becomes unstable through 
a subcritical pitchfork bifurcation, due to its collision 
with two unstable asymmetric branches. The VA pre- 
dicts this critical point at E ^ 3.92, in good agreement 
with the numerical findings. Further, at another criti- 
cal value, E = 3.166 (the corresponding VA prediction is 
E ^ 3.128) the unstable asymmetric branches turn back 



as stable ones, through a saddle-node bifurcation. Be- 
tween the two critical points, both the symmetric branch 
and the outer asymmetric one are stable, hence there ex- 
ists a region of bistability. To additionally demonstrate 
the accuracy of the VA, Fig. [4] presents comparison of the 
solution profiles dit E = 3.4, together with the spectral 
plane, (Re(A), Im(A)), for the corresponding eigenvalues, 
A = Re(A) +iIm(A). Recall that in Hamiltonian systems, 
such as the one considered here, if A is an eigenvalue, so 
are also —A, A"^ and — A"*", hence, if an eigenvalue with a 
nonzero real part exists, then the system will be linearly 
unstable. 

For those ID solutions that are linearly unstable due 
to a real eigenvalue, such as the unstable asymmetric so- 
lutions for 3.166 < E < 3.82 and the symmetric ones for 
E > 3.82, we have examined their evolution in Fig. [5] by 
means of direct simulations of Eq. ([1]). In the case of the 
unstable asymmetric branch, the result of the evolution 
depends on the nature of the corresponding perturbation, 
due to the presence of the bistability in the corresponding 
parameter range: the solution ends up oscillating around 
either the stable asymmetric solution, or the stable sym- 
metric one. The evolution of the unstable symmetric 
branch naturally results in oscillations around the stable 
asymmetric profile, which represents the ground state in 
that case. 

An important observation in comparing Figs. [1] and 
[2] is that the bifurcation found in the anti-continuum 
limit shown in Fig. [1] is definitely supercritical, unlike the 
weakly subcritical one in Fig. [2l This indicates that the 
character of the bifurcation changes from subcritical to 
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FIG. 8: Time-evolution plots for a perturbation of the unstable stationary state from Fig. [6] The top-row perturbation pushes 
the solution towards the upper stable part of the curve, while the bottom-row perturbation pushes it to the symmetric branch. 
The left column shows the three-dimensional space-time evolution of isodensity contours of the U (respective top subpanels) 
and V (respective bottom subpanels) solutions, while the right column shows the evolution of the asymmetry measure, r, from 
an unstable steady state towards a stable one (both are denoted by dashed lines). 



supercritical pitchfork with the increase of discreteness, 
i.e., decrease of e. This transition should eliminate the 
unstable asymmetric branches. In accordance with this 
expectation, we have found that the unstable asymmetric 
solutions exist only for e > 0.35, in the ID case. 

We now turn to the results for the 2D model collected 
in Fig. [3l Here, the results are even more interesting, 
for a number of reasons. On the one hand, there is no 
continuum analog to the bifurcation diagram, as the re- 
spective 2D continuum solutions are always unstable to 
collapse. Discreteness is well-known to arrest collapse 
[22! , [27! , generating branches of potentially stable local- 
ized solutions for sufficiently small values of the inter- 
site coupling constant (at a given chemical potential), or 
for sufficiently large chemical potential (at a given inter- 
site coupling), in one-component models [22]. Further- 
more, in the 2D case, for a given value of the norm, there 
are two coexisting symmetric solutions, one (taller and 
narrower) with a larger chemical potential, which is sta- 
ble, and one (shorter and wider) with a smaller chemi- 
cal potential, which is unstable. As Fig. [3] implies, the 
symmetry-breaking weakly subcritical pitchfork bifurca- 
tion typically occurs from the stable branch of the sym- 
metric solution, the corresponding critical point in Fig. [3] 
being E ^ 1.45. Similarly to the ID case, there is also a 
saddle-node bifurcation between the unstable and stable 
asymmetric branches, which occurs at E ^ 1.411 and is 



responsible for the turning point. 

In the 2D model too, the VA accurately captures the 
trends of the numerical results, even though the less ac- 
curate nature of the 2D ansatz prevents a quantitative 
matching of the resulting bifurcation diagrams. Detailed 
profiles of the numerical solutions and their VA-predicted 
counterparts are shown in Figs. [6] and [71 for the asymmet- 
ric and symmetric branches respectively, dit E = 1.435. 
Two additional remarks are in order here. Firstly, simi- 
lar to the ID case, the 2D bifurcation changes character 
from weakly subcritical (as observed in Fig. [3]) to super- 
critical (as seen in the anti-continuum limit of e = 0) at 
e ~ 0.19. On the other hand, as e further increases, the 
stable part of the symmetric branch (in Fig. [3]) shrinks 
and eventually disappears at e > 0.29. 

Finally, we have again examined the dynamics of lin- 
early unstable solutions, upon appropriate perturbations, 
through direct simulations, see Figs. [8] and [9l In the for- 
mer figure, we have explored how the bist ability, which 
is shown for a range of E in Fig. O "kicks" the unstable 
asymmetric solution either in the direction of its stable 
asymmetric counterpart, or towards the stable symmet- 
ric solution [it is relevant to stress here, in connection 
with Fig. O that, while the symmetric solution has a 
norm threshold at £^ ~ 1.425, the stable asymmetric so- 
lution can, in principle, be found at lower norms than the 
symmetric one, thus allowing the system to effectively de- 
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FIG. 9: Similar to the previous figure: the time- evolution plots for perturbations of two unstable symmetric steady states. The 
top row is for a solution with ^ = 1.5, and the bottom row is for the solution from Fig. [T] In the top row, the perturbation 
pushes the solution towards the upper asymmetric branch, while in bottom row the perturbation initiates the evolution of the 
solution towards the stable symmetric state. 



crease its "excitation threshold" [28^]. In the latter figure, 
we consider the unstable symmetric branch, both when 
a stable symmetric branch does not exist, in which case 
the solution becomes asymmetric (top row), and when 
a stable symmetric branch does exist, in which case the 
system evolves towards that solution (bottom row). 

IV. CONCLUSIONS 

In this work, we have introduced the model based on 
two linearly coupled lattices with the cubic nonlinearity, 
and investigated its dynamical properties in detail. We 
have demonstrated that, in a number of respects, the 
discrete system emulates its continuum ID counterpart, 
analyzed earlier in Refs. [11,1135 118]. On the other hand, 
the lattice model gives rise to novel features. Even in 
the ID setting, varying the strength of the lattice inter- 
site coupling may be used to switch the character of the 
bifurcation from subcritical to supercritical pitchfork, as 
the coupling gets weaker, and the anti-continuum limit is 
approached. In the more interesting 2D setting, our work 
is the first manifestation, to the best of our knowledge, of 
the existence of such a bifurcation diagram, since in the 



continuum limit the symmetry-breaking models with the 
self-focusing nonlinearity are irrelevant due to the col- 
lapse instability. Furthermore, the discreteness induces 
the presence of both stable and unstable branches of sym- 
metric solutions, thus enriching the bifurcation diagram. 
In the 2D case, not only is it possible for weaker lattice 
coupling to turn the bifurcation from subcritical to super- 
critical, but it is also possible for the lattice (when the 
bifurcation is subcritical) to possess a lower excitation 
threshold for asymmetric states than for the symmetric 
ones. All of these features demonstrate critical modifi- 
cations that the discreteness imposes on the well-known 
symmetry-breaking picture in continuum models. 

It might be quite interesting to examine similar fea- 
tures in 3D and compare the results with their 2D and 
ID counterparts. Of perhaps even more physical interest, 
especially in terms of coupled hyperfine states in BECs, 
would be to add nonlinear coupling between the lattice 
components, of the cross-phase-modulation type, to the 
linear coupling considered here. It would be particularly 
interesting to examine how the symmetry-breaking phe- 
nomenology is affected by the gradual increase of such 
a coupling. This study is currently in progress and the 
results will be reported elsewhere. 
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